/*******************************************************************************

                                WEEDFC.C

  Title:        Speden weed proedure for fcalc files. 
  Author:       Hanna Szoke
  Date:         12/30/2004 
  Function:     Weed out near-duplicate structure factors.
		It is assumed that the file has previously been sorted,
		at least to the extent that the h values are clumped,
		using, e.g.  'sort -n +1 -2 [in_filename] >[out_filename]',
		or even just 'sort [in_filename] >[out_filename]'.

*******************************************************************************/

#include	"util.h"
#include	"cellparams.h"
#include	"dims.h"
#include	"sp.h"

float	fcround = 0.1 ;
char	weed_filename[MAXSTRING], dups_filename[MAXSTRING] ;
FILE	*fp_weed, *fp_dups ;

void
weedfc_main(int argc, char *argv[]) {

     int	HDFdatalength(char *, char *) ;
     int	readFcList(char *) ;
     int	read_HDFfcalc(char *) ;
     void	alloc_fc(int) ;
     int	weedit_fc(int) ;
     void	writeFclist(char *) ;
     void	write_HDFfcalc(char *) ;

     char	infilename[MAXSTRING] ;
     char	outfilename[MAXSTRING] ;
     int	NoDup ;

     if (argc > optind+5) 	/* too many arguments */
	  ballpark(caller) ; 
     else if (argc > optind+3) {
          sprintf(infilename, "%s", argv[optind+2]) ;
          sprintf(outfilename, "%s", argv[optind+3]) ;
          if (argc == optind+5)
               fcround  = (float) atof(argv[optind+4]) ;
     } 
     else {
	  prompt("\nPlease enter 2 filenames - input and output: ") ;
          sscanf(terminp, "%s  %s", infilename, outfilename) ;
     }

     sprintf(weed_filename, "%s.weed", infilename) ;
     sprintf(dups_filename, "%s.dups", infilename) ;

     if ((fp_weed = fopen(weed_filename,"w")) == NULL) {
          sprintf(message, "Cannot open %s", weed_filename);
          Fatal(message) ;
     }
     if ((fp_dups = fopen(dups_filename,"w")) == NULL) {
          sprintf(message, "Cannot open %s", dups_filename);
          Fatal(message) ;
     }

     hello(caller) ;

     readBasicInput() ;
     if (HDFIN)
          Nfc = HDFdatalength(infilename, "/fcalc/Nfc") ;
     else
          Nfc = filelength(infilename) ;

     strcpy(message, "weedfc") ;

     alloc_fc(Nfc) ;
     printTwice("Reading input...") ;
     if (HDFIN)
          Nfc = read_HDFfcalc(infilename) ;
     else
          Nfc = readFcList(infilename) ;
     printTwice("Weeding...") ;
     NoDup = weedit_fc(Nfc) ;
     Nfc = NoDup ;
     printTwice("Writing output.") ;
     if (HDFOUT)
          write_HDFfcalc(outfilename) ;
     else
          writeFclist(outfilename) ;
}
int weedit_fc(N)
{

	/****************************************************************
	Check each entry against all others.  Is it a (near) duplicate 
	in (h,k,l)?  (Leaving duplicates means weighting them 
	disproportionately.)  If so, replace all such by the mean.
	Report duplicates and statistics of weeding.
	****************************************************************/

     real	cleanPhase(real *) ;
     float	roundit(float *, float) ;
     int        n, nn ;
     int	ndups, totdups = 0 ;
     float	h, k, l ;
     float	hh, kk, ll ;
     real	nextre, nextim ;
     real 	avere, aveim, relstd ;
     real	omega, amplitude, phase ;

     int	*dups ;
     char	*weeds ;

     dups = (int *) e_malloc(N*sizeof(int), "weed_duplicates") ;
     weeds = (char *) e_malloc(N*sizeof(char), "weed_duplicates") ;

     for (n = 0; n < N; n++) {
	  *(dups+n) = 1 ;
          *(weeds+n) = 0 ;
     }

		/* Round all indices to multiples of 'fcround'	*/

     for (n = 1; n < N; n++) {
        *(Fc_Hlist+n) = roundit(Fc_Hlist+n, fcround) ;
        *(Fc_Klist+n) = roundit(Fc_Klist+n, fcround) ;
        *(Fc_Llist+n) = roundit(Fc_Llist+n, fcround) ;
     }

	/* Identify entries that have duplicates ***/

     for (n = 0; n < N; n++) {

        if (*(weeds+n) == 0) {
           h = *(Fc_Hlist + n) ;
           k = *(Fc_Klist + n) ;
           l = *(Fc_Llist + n) ;

           for (nn = n+1; nn < N; nn++) {
   
              hh = *(Fc_Hlist + nn) ;
              kk = *(Fc_Klist + nn) ;
              ll = *(Fc_Llist + nn) ;

              if ((h == hh) && (k == kk) && (l == ll)) {
                 *(dups+n) += 1 ;
                 *(weeds+nn) = 1 ;
                 totdups += 1 ;
              }
              else if (hh != h) 
                 break ;

           }
        }
     }

	/* determine what should replace the duplicates ***/

     for (n = 0; n < N; n++) {

        if (*(dups+n) > 1) {
           ndups = *(dups+n) ;
           h = *(Fc_Hlist + n) ;
           k = *(Fc_Klist + n) ;
           l = *(Fc_Llist + n) ;

           avere = *(Fc_relist + n) ;
           aveim = *(Fc_imlist + n) ;
           relstd = avere*avere + aveim*aveim ;

           amplitude = sqrt(avere*avere + aveim*aveim) ; 
	   omega = RTOD * atan2(aveim, avere) ;
           phase = cleanPhase(&omega) ;
	   fprintf(fp_dups, "INDEX   %10.2f  %10.2f  %10.2f  FCALC   %g   %g\n", 
	          h, k, l, amplitude, phase) ;

           for (nn = n+1; nn < N; nn++) {

              hh = *(Fc_Hlist + nn) ;
              kk = *(Fc_Klist + nn) ;
              ll = *(Fc_Llist + nn) ;

              if ((h == hh) && (k == kk) && (l == ll)) {
                 nextre = *(Fc_relist + nn) ;
                 nextim = *(Fc_imlist + nn) ;
                 avere += nextre;
                 aveim += nextim;
                 relstd += nextre*nextre + nextim*nextim ;

                 amplitude = sqrt(nextre*nextre + nextim*nextim) ; 
	         omega = RTOD * atan2(nextim, nextre) ;
                 phase = cleanPhase(&omega) ;
	         fprintf(fp_dups, 
                         "INDEX   %10.2f  %10.2f  %10.2f  FCALC   %g   %g\n", 
	  	         hh, kk, ll, amplitude, phase) ;
              }
              else if (hh != h) 
                 break ;
           }
           avere = avere / ndups ;
           aveim = aveim / ndups ;
           if ((avere != 0) || (aveim != 0)){
                relstd -= ndups * (avere*avere + aveim*aveim) ;
                relstd = sqrt(fabs(relstd) / (ndups-1))/ 
                         sqrt(avere*avere + aveim*aveim) ;
          }
          else
                relstd = 0 ;

          fprintf(fp_weed, 
             "%d inputs near (%g,%g,%g) replaced by %g, %g rel std: %g \n", 
		ndups, h, k, l, avere, aveim, relstd) ;

           *(Fc_relist + n) = avere ;
           *(Fc_imlist + n) = aveim ;
         }
      }
      fclose(fp_weed) ;
      fclose(fp_dups) ;

     /*** Finally - do the weeding ***/

     n = 0 ;
     for (nn = 0; nn < N; nn++) {
         if (*(weeds+nn) == 0) {
            *(Fc_Hlist + n) = *(Fc_Hlist + nn) ;
            *(Fc_Klist + n) = *(Fc_Klist + nn) ;
            *(Fc_Llist + n) = *(Fc_Llist + nn) ;

            *(Fc_relist + n) = *(Fc_relist + nn) ;
            *(Fc_imlist + n) = *(Fc_imlist + nn) ;
            n++ ;
        }
     }

     sprintf(message, "A total of %d inputs have been removed", totdups) ;
     printTwice(message) ;
     fprintf(stdout, "For details, see .weed and .dups files\n\n") ;

     e_free(N*sizeof(int), dups) ;
     e_free(N*sizeof(char), weeds) ;

     return(n) ;
}

